#include <cstdio>
#include <cstdlib>
#include <cstring>
#include "parse.h"
#include <cuda_runtime.h>
#include "texArrayHandle.h"
#include <vector>
#include "molecule.h"
#include "energy_ruler.h"

using std::vector;

extern const int IDX_E;
extern const int IDX_D;
extern const int IDX_VDW;
extern const int IDX_SOL;
extern const int IDX_EPS;
extern const int IDX_R_EPS;

extern   void init_texture(tex3DHandle<float>&, int, bool);
extern   void init_texture(tex1DHandle<float>&, int);
//extern "C" void launch_kernel(float *, float, float, float);
//extern "C" void launch_kernel(int);
extern "C" void launch_kernel(int, float, float, float, float, float, float ,float, float *, short *);
//extern "C" void launch_kernel(int, short *);
extern "C" void ligandCopy(ligand &);
extern "C" void mapSet(float &, float *);

int main(int argc, char * argv[])
{
	// SCAN DPF
	if (argc!=2)
	{
		printf("Please specify a valid DPF file\n");
		return 1;
	}
	FILE * DPF = fopen(argv[1],"r");
	if (DPF==NULL)
	{
		printf("Please specify a valid DPF file\n");
		return 1;
	}

	char line[1000], temp[1000], var[1000];
	char lig_name[1000], flex_name[1000]={};
	int dim[3],mapCount;
	char mapName[1000][50];
	int skip[50];
	float spacing=0;
	float max_coord[3];
	float min_coord[3];
	float center[3];
	
	while (fgets(line,1000,DPF)!=NULL)
	{
		int token=parse_dpf(line);
		switch (token)
		{
			case DPF_fld:
			{
				sscanf(line,"%s %s", temp, var);
				FILE * FLD = fopen(var,"r");
				FILE * XYZ;

				// LOAD MAPS INFO
				/*
				   1. scan maps.fld for dims, map names and line skipping
			   	   2. scan maps.xyz for the boundary
				*/
				while (fgets(line,1000,FLD)!=NULL)
				{
					if (line[0]=='#') continue;
					if (strncmp(line,"dim1",4)==0)
					{
						sscanf(line,"%4s=%d",temp,&dim[0]);
					}
					else if (strncmp(line,"dim2",4)==0)
					{
						sscanf(line,"%4s=%d",temp,&dim[1]);
					}
					else if (strncmp(line,"dim3",4)==0)
					{
						sscanf(line,"%4s=%d",temp,&dim[2]);
					}
					else if (strncmp(line,"veclen",6)==0)
					{
						sscanf(line,"%6s=%d",temp,&mapCount);
					}
					else if (strncmp(line,"coord 1",7)==0)
					{
						sscanf(line,"%s %s %4s=%s ",temp,temp,temp,var);
					}
					else if (strncmp(line,"variable",8)==0)
					{
						int index;
						sscanf(line,"%s %d",temp,&index);
						sscanf(line,"%s %s %4s=%s %s %4s=%d",temp,temp,temp,mapName[index-1],temp,temp,&skip[index-1]);
					}
				}
				
				XYZ=fopen(var,"r");
				for (int i=0;i<3;i++)
				{
					fscanf(XYZ,"%f %f",&min_coord[i],&max_coord[i]);
					center[i]=((min_coord[i]+max_coord[i])/2.0);
				}

				fclose(XYZ);


				fclose(FLD);
				break;
			}
			case DPF_map:
				break;
			case DPF_move:
				sscanf(line,"%*s %s\n",lig_name);
				break;
			case DPF_flex:
				sscanf(line,"%*s %s\n",flex_name);
				break;
			default:
				break;
		}
	}

	fclose(DPF);
	// END SCAN DPF

	printf("Grid configurations:\n");
	printf("center\t[%f %f %f]\n",center[0],center[1],center[2]);
	printf("min\t[%f %f %f]\n",min_coord[0],min_coord[1],min_coord[2]);
	printf("max\t[%f %f %f]\n",max_coord[0],max_coord[1],max_coord[2]);
	printf("dims\t[%d %d %d]\n",dim[0],dim[1],dim[2]);
	spacing=(max_coord[0]-min_coord[0])/((float)dim[0]-1);
	printf("spacing\t%.4f\n\n",spacing);

	mapSet(spacing,min_coord);

	// LOADING MAPS
	FILE * MAP;
//	float maps[mapCount][dim[2]][dim[1]][dim[0]];	// stack may not be able to hold this much
	float * maps=(float * )malloc(sizeof(float)*mapCount*dim[2]*dim[1]*dim[0]);
/*	float **** maps=(float ****)malloc(mapCount*sizeof(float ***));
	for (int i=0;i<mapCount;i++)
	{
		maps[i]=(float ***)malloc(dim[2]*sizeof(float **));
		for (int j=0;j<dim[2];j++)
		{
			maps[i][j]=(float **)malloc(dim[1]*sizeof(float *));
			for (int k=0;k<dim[1];k++)
			{
				maps[i][j][k]=(float *)malloc(dim[0]*sizeof(float));
			}
		}
	}
*/		
	tex3DHandle<float> tex[mapCount];
//	vector< tex3DHandle<float> > tex(mapCount);

	for (int i=0;i<mapCount;i++)
	{
//		if (i!=7) continue;
		printf("Map #%d: %s\n",i,mapName[i]);
		MAP=fopen(mapName[i],"r");
		for (int j=0;j<skip[i];j++)
			fgets(line,1000,MAP);
		for (int z=0;z<dim[2];z++)
			for (int y=0;y<dim[1];y++)
				for (int x=0;x<dim[0];x++)
				{
					fscanf(MAP,"%f",&maps[i*dim[2]*dim[1]*dim[0]+z*dim[1]*dim[0]+y*dim[0]+x]);
				}
		fclose(MAP);

		tex3DHandle<float> temp(dim[0],dim[1],dim[2]);

//		tex.push_back(temp);
		tex[i]=temp;
//		tex[i].print();
		tex[i].copyData((void*)(maps+i*dim[2]*dim[1]*dim[0]));
		init_texture(tex[i],getIDX(mapName[i]),true);
	}

	printf("\n");

	// LOADING ligand
	//
	
	int atomCount=0;
	int torCount=0;
/*
	if (strlen(flex_name)!=0)
	{
*/
		/*
		ligand flex(flex_name);
		printf("flex\n");
		flex.print_tor();
		*/

//		ligand moving(lig_name,flex_name);
		ligand moving(lig_name);
		printf("total\n");
		moving.print_tor();
		ligandCopy(moving);
		atomCount=moving.atomCount;
		torCount=moving.torCount;
/*
		for (int i=0;i<atomCount;i++)
		{
			printf("%8d  |",i);
			for (int j=0;j<atomCount;j++)
			{
				if (moving.nonbonded[i][j]==0)
					printf("_|");
				else printf("X|");
			}
			printf("\n");
		}
*/
		energyRuler ruler;

/*	
		for (int i=0;i<1000;i++)
		{
			//printf("vdw [5][5][%3d] = %8.3f\n",i,ruler.vdw_table[5*27*1000+5*1000+i]);
			printf("r_epsilon [%3d] = %8.3f\n",i,ruler.r_epsilon[i]);
			printf("desol [%3d] = %8.3f\n",i,ruler.sol_table[i]);
			printf("epsilon [%3d] = %8.3f\n",i,ruler.epsilon[i]);
		}
*/		
//		ruler.print();
//
		tex3DHandle<float> vdw_map(1000,27,27);
		
		vdw_map.copyData((void*)(ruler.vdw_table));
		init_texture(vdw_map,IDX_VDW,false);

		tex1DHandle<float> desol(1000);
		desol.copyData((void*)(ruler.sol_table));
		init_texture(desol,IDX_SOL);

		tex1DHandle<float> epsilon(1000);
		epsilon.copyData((void*)(ruler.epsilon));
		init_texture(epsilon,IDX_EPS);

		tex1DHandle<float> r_eps(1000);
		r_eps.copyData((void*)(ruler.r_epsilon));
		init_texture(r_eps,IDX_R_EPS);
		/*
	}
	else
	{
		ligand move(lig_name);
		move.print_tor();
		ligandCopy(move);
		atomCount=move.atomCount;
		torCount=move.torCount;
	}
	*/
/*
	float coord[]={34.602669,29.242666,40.922668};
	
	float *d_out;
	cudaMalloc(&d_out,sizeof(float)*10);
	float *h_out;
	h_out=(float*)malloc(sizeof(float)*10);
	
	launch_kernel(d_out,coord[0],coord[1],coord[2]);

	cudaMemcpy(h_out,d_out,sizeof(float)*10,cudaMemcpyDeviceToHost);
	printf("\nGetTexture\n");
	for (int i=0;i<8;i++)
	{
		printf("%.5f\n",h_out[i]);
	}

	printf("Int\t%.5f\n",h_out[8]);
	printf("Tex\t%.5f\n",h_out[9]);
*/
	float angle[]={167.35,-44.14,102.58,-162.99,144.17,59.53,133.93,-49.18,139.08,23.82,-0.88,70.10,-72.69};
//	float angle[]={0,0,0,0,0,0,0,0,0,0,0,0,0};
	float h_theta_half[torCount];
	float *theta_half;
	cudaMalloc(&theta_half,sizeof(float)*torCount);

	//0.0174532925=(pi/180)
	//0.00872664626 = (pi/180)/2
	for (int i=0;i<torCount;i++)
	{
		h_theta_half[i]=angle[i]*0.00872664626f;
		printf("(%6.3f)",h_theta_half[i]);
	}

	cudaMemcpy(theta_half,h_theta_half,sizeof(float)*torCount,cudaMemcpyHostToDevice);
	printf("\n");	
//	printf("\n[%.4f %.4f %.4f]\n",move.cent_x,move.cent_y,move.cent_z);

//	launch_kernel(atomCount,0.368900,-0.214800,-4.986500,0,0,0,1,theta_half);
//
	short * nonBond;
	printf("atomCount %d\n",atomCount);
//	size_t pitch;
//	cudaMallocPitch((void **)&nonBond,&pitch,atomCount*sizeof(short),atomCount);
	cudaMalloc((void**)&nonBond,sizeof(short)*atomCount*atomCount);

//	printf("pitch \t %ld\n",pitch);

	short tempMatrix[atomCount*atomCount];
	for (int i=0;i<atomCount;++i)
		for (int j=0;j<atomCount;++j)
		{
			tempMatrix[i*atomCount+j]=moving.nonbonded[i][j];
		}

//	cudaMemcpy2D(nonBond,pitch,tempMatrix,atomCount,atomCount,atomCount,cudaMemcpyHostToDevice);
	cudaMemcpy(nonBond,tempMatrix,sizeof(short)*atomCount*atomCount,cudaMemcpyHostToDevice);

	launch_kernel(atomCount,2.543500,6.025500,-7.966300,0.529539,-0.490877,-0.339952,0.602545,theta_half,nonBond);
	//launch_kernel(atomCount,nonBond);
	return 0;
}
